;+
; NAME:
;       OMNI_COMPUTE_PHYSICS
;
; PURPOSE:
;       Compute physical properties and propogate uncertainties in a
;       Monte Carlo fashion for SURVEY catalog objects.
;
; CATEGORY:
;       distance-omnibus Mass Routine
;
; CALLING SEQUENCE:
;       OMNI_COMPUTE_PHYSICS
;
; INPUTS:
;       NONE
;
; OPTIONAL INPUTS:
;       NONE
;
; KEYWORD PARAMETERS:
;       NONE
;
; OUTPUTS:
;       NONE [Saves physics to an IDL save file.]
;
; OPTIONAL OUTPUTS:
;       NONE
;
; COMMON BLOCKS:
;       OMNI_CONFIG -- The set of configuration structures, read in
;                      from the config files in conffiles/
;       MASS_BLOCK  -- Set of variables read in and computed in
;                      OMNI_MASS_SETUP which are used in other
;                      mass-related routines.
;
; MODIFICATION HISTORY:
;
;       Created:  02/04/14, TPEB -- Initial version.
;       Modified: 02/27/14, TPEB -- Attempts at speed-up.
;
;-

PRO OMNI_COMPUTE_PHYSICS, CONFFILE=cfile, SAMPLETD=sampletd, VERBOSE=verbose, $
                          SILENT=silent
  
  COMPILE_OPT IDL2, LOGICAL_PREDICATE, HIDDEN
  
  COMMON OMNI_CONFIG, conf, mw, local, dpdfs, ancil, fmt, conffile
  COMMON MASS_BLOCK, pvec, s, iflux, ieflux, m0, hnuk, kappa, nkp
  
  ;; Parse keywords
  verbose  = KEYWORD_SET(verbose) && ~KEYWORD_SET(silent)
  silent   = KEYWORD_SET(silent)  && ~KEYWORD_SET(verbose)
  
  ;; Read in survey-info, galactic-params, & dpdf-params config files
  conf     = omni_load_conf(cfile)
  mw       = omni_read_conffile('./conffiles/galactic_params.conf')
  dpdfs    = omni_read_conffile('./conffiles/dpdf_params.conf')
  mconf    = omni_read_conffile('./conffiles/mass_deriv.conf')
  d        = dindgen(dpdfs.nbins)*dpdfs.binsize + dpdfs.binstart
  sampletd = n_elements(sampletd) ? sampletd GE 1 : mconf.sampletd
  
  ;; Load up the various COMMON block variables using
  ;;   configuration-file instructions.
  status = OMNI_MASS_SETUP(CONFFILE=cfile, VERBOSE=verbose, SILENT=silent, $
                          TYPSTR='physical property')
  IF ~status THEN RETURN
    
  ;; Read in the SURVEY catalog 
  restore,'./local/'+conf.survey+'_map_locations.sav',VERBOSE=verbose
  
  ;; Limit SURVEY structure to the items saved in the above line
  surv_ind = WHERE_ARRAY(s.cnum, survey.cnum, n_survey)
  IF n_survey NE n_elements(s) THEN message,'Rum now, Think later!!!'
  survey = survey[surv_ind]
  
  
  ;;=========================================================================
  ;; Monte Carlo the mass based on the DPDF, Gaussian distribution of
  ;;   flux density, and (if desired) lognormal distribution of
  ;;   temperature.
  
  ndloop = 1d7                  ; TEN MILLION iterations!
  sampletd=1b                   ; Force temperature sampling
  
  min_arad = min(s.rad, /NaN)   ; Define minimum angluar radius for checking
  min_arad = 5.

  ;;===============&&&&&&&&&&&&&&&===============
  ;; Create a structure to contain the N = large random-sampled values
  ;;   of each of distance, flux density, and tempertature.
  obj = {dist: dblarr(ndloop), $
         sint: dblarr(ndloop), $
         temp: dblarr(ndloop) + mw.td, $
         mass: dblarr(ndloop), $
         m20k: dblarr(ndloop), $
         prad: dblarr(ndloop), $
         nden: dblarr(ndloop), $
         n20k: dblarr(ndloop)}
  
  ;; Create a structure to contain the output ML values and
  ;;    standardized PDFs.
  phys = replicate( {cnum: 0L,$
                     glon: 0.d,$
                     glat: 0.d,$
                     rgal: 0.d,$
                     zgal: 0.d,$
                     mass: fltarr(3), $    ; Log10 version
                     m20k: fltarr(3), $    ; Fix T_d = 20K
                     prad: fltarr(3), $    ; Log10 version
                     nden: fltarr(3), $    ; Log10 version
                     n20k: fltarr(3), $    ; Fix T_d = 20K
                     mpdf: fltarr(mconf.nmbin), $
                     m20f: fltarr(mconf.nmbin), $
                     rpdf: fltarr(mconf.nrbin), $
                     npdf: fltarr(mconf.nnbin)}, nkp)
  
  ;; Compute the abscissae for the physical quantity PDFs
  m_absc = findgen(mconf.nmbin) * mconf.mbinsize + mconf.mbinstart
  r_absc = findgen(mconf.nrbin) * mconf.rbinsize + mconf.rbinstart
  n_absc = findgen(mconf.nnbin) * mconf.nbinsize + mconf.nbinstart
  ;;===============&&&&&&&&&&&&&&&===============
  
  
  IF sampletd THEN BEGIN
     ;; If sampletd is set, set up what is needed for providing
     ;;   a random temperature distribution.
     ;;   Temperature distribution is LOG-NORMAL!!!
     ;; 0K - 100K base array in 0.02K steps
     tdarr = dindgen(5001)/50. 
     lnsig = ( alog(mconf.tdmu + 0.5*mconf.tdfwhm) - $
               alog(mconf.tdmu - 0.5*mconf.tdfwhm) ) / 2.355
     tdpdf = gauss_1(alog(tdarr),[1.,alog(mconf.tdmu),lnsig])
  ENDIF
  
  message,'Computing physics... please hold for the next available customer '+$
          'service representative.',/inf
  print,''
  
  openw,lun,'./masses/mtable.tex',/get_lun
  openw,lun2,'./masses/mtable2.tex',/get_lun
  
  ;; Create arrays for analysis of proper PDF abscissae ranges.
  mmax = dblarr(nkp)
  rmax = dblarr(nkp)
  nmax = dblarr(nkp)
  mmin = dblarr(nkp)
  rmin = dblarr(nkp)
  nmin = dblarr(nkp)
  
  ;; Load the outout structure elements
  phys.cnum = pvec.cnum
  phys.glon = pvec.glon
  phys.glat = pvec.glat
  omni_lbd2rz,phys.glon,phys.glat,pvec.stat.duse[0],Rgal,Z
  phys.rgal = Rgal / 1.d3
  phys.zgal = Z
  
  ;; Time this POS routine.
  start_t = systime(1)
  
  ;; Loop through objects!
  FOR ii=0, nkp-1 DO BEGIN
     
     message,'Object '+string(ii+1,format="(I0)")+' of '+$
             string(nkp,format="(I0)"),/inf
     
     ;; Load up distances based on this posterior DPDF
     obj.dist = omni_sample_dpdf(pvec[ii].post,d,ndloop)
     
     ;; Compute the flux density PDF, then load up fluxes!
     sran = s[ii].(iflux) + [-1.d,1.d]*s[ii].(ieflux)*10   ; 10-sigma
     dels = sran[1] - sran[0]                              ; Width of above
     sarr = dindgen(1001)/1000. * dels + sran[0]           ; 1000-element arr
     spdf = gauss_1(sarr,[1.,s[ii].(iflux),s[ii].(ieflux)]) ; PDF
     sind = where(sarr LT survey[ii].noise, nsi)
     IF nsi NE 0 THEN spdf[sind] = 0. ; Set probability = 0 for S < noise
     obj.sint = omni_sample_dpdf(spdf,sarr,ndloop)
     
     ;; If desired, load random temperatures
     IF sampletd THEN obj.temp = omni_sample_dpdf(tdpdf,tdarr,ndloop)
     
     ;; Compute temperature correction; == 1 if sampletd = 0b.
     tempcorr = (exp(hnuk/obj.temp)-1.d) / $
                (exp(hnuk/mw.td)-1.d)
     
     ;;###################################################################
     ;; Compute the physics!!!
     
     ;; Mass in M_sun
     obj.mass = m0 * obj.sint * conf.fluxcor * tempcorr * $
                obj.dist * obj.dist / 1.d6
     obj.m20k = m0 * obj.sint * conf.fluxcor * $
                obj.dist * obj.dist / 1.d6

     ;; Physical radius in pc
     IF ~finite(s[ii].rad) THEN BEGIN
        s[ii].rad = min_arad    ; Check angular radius
        thetar = '\nodata'
        unresolved = 1b
     ENDIF ELSE BEGIN
        thetar = string(s[ii].rad,format="(F0.1)")
        unresolved = 0b
     ENDELSE
     
     obj.prad = (s[ii].rad / 206265.) * obj.dist
     
     ;; Number density in cm^{-3}
     obj.nden = obj.mass * !pc.msun / (!pc.mu * !pc.mh) / $
                ( !dpi * (obj.prad * obj.prad * obj.prad * $
                          !pc.pc * !pc.pc * !pc.pc * 1.d6) )
     obj.n20k = obj.m20k * !pc.msun / (!pc.mu * !pc.mh) / $
                ( !dpi * (obj.prad * obj.prad * obj.prad * $
                          !pc.pc * !pc.pc * !pc.pc * 1.d6) )
     
     ;; Convert arrays to LOG10 for ease of computing things...
     obj.mass    = alog10(obj.mass)
     obj.m20k    = alog10(obj.m20k)
     obj.prad    = alog10(obj.prad)
     obj.nden    = alog10(obj.nden)
     obj.n20k    = alog10(obj.n20k)
     
     
     ;;###################################################################
     ;; Statistics, baby!
     
     ;; ===== Mass =====
     plothist,obj.mass,bin=mconf.mbinsize,marr,mhis,/noplot
     phys[ii].mass = omni_compute_ebars(mhis,marr,/logcut,NJ=nj)
     IF NJ GE 2 THEN message,'Fire 1',/inf
     mmax[ii] = max(marr)
     mmin[ii] = min(marr)
     
     marr = [(findgen(10)-10.)*mconf.mbinsize + mmin[ii],marr,$ ; Pad the raw 
             (findgen(10)+1)*mconf.mbinsize + mmax[ii]]         ; histogram for
     mhis = [fltarr(10),mhis,fltarr(10)]                        ; interpolation
     
     phys[ii].mpdf = interpol(mhis,marr,m_absc) ; Interpolate onto common scale
     phys[ii].mpdf /= total(phys[ii].mpdf) ; Normalize to unit integral prob.
     
     
     plothist,obj.m20k,bin=mconf.mbinsize,marr,mhis,/noplot
     phys[ii].m20k = omni_compute_ebars(mhis,marr,/logcut,NJ=nj)
     IF NJ GE 2 THEN message,'Fire 2',/inf
     mmax[ii] = max(marr)
     mmin[ii] = min(marr)
     
     marr = [(findgen(10)-10.)*mconf.mbinsize + mmin[ii],marr,$ ; Pad the raw 
             (findgen(10)+1)*mconf.mbinsize + mmax[ii]]         ; histogram for
     mhis = [fltarr(10),mhis,fltarr(10)]                        ; interpolation
     
     phys[ii].m20f = interpol(mhis,marr,m_absc) ; Interpolate onto common scale
     phys[ii].m20f /= total(phys[ii].m20f) ; Normalize to unit integral prob.
     
     
     ;; ===== Physical Radius =====
     plothist,obj.prad,bin=mconf.rbinsize,parr,phis,/noplot
     phys[ii].prad = omni_compute_ebars(phis,parr,/logcut,NJ=nj)
     rmax[ii] = max(parr)
     rmin[ii] = min(parr)
     
     parr = [(findgen(10)-10.)*mconf.rbinsize + rmin[ii],parr,$ ; Pad the raw 
             (findgen(10)+1)*mconf.rbinsize + rmax[ii]]         ; histogram for
     phis = [fltarr(10),phis,fltarr(10)]                        ; interpolation
     
     phys[ii].rpdf = interpol(phis,parr,r_absc) ; Interpolate onto common scale
     phys[ii].rpdf /= total(phys[ii].rpdf) ; Normalize to unit integral prob.
     cgPlot,r_absc,phys[ii].rpdf,/xst
     IF NJ GE 2 THEN message,'Fire 3',/inf
     
     
     ;; ===== Number Density =====
     plothist,obj.nden,bin=median(obj.nden)/1.d3,narr,nhis,/noplot
     phys[ii].nden = omni_compute_ebars(nhis,narr,/logcut,NJ=nj)
     nmax[ii] = max(narr)
     nmin[ii] = min(narr)
     
     narr = [(findgen(10)-10.)*mconf.nbinsize + nmin[ii],narr,$ ; Pad the raw 
             (findgen(10)+1)*mconf.nbinsize + nmax[ii]]         ; histogram for
     nhis = [fltarr(10),nhis,fltarr(10)]                        ; interpolation
     
     phys[ii].npdf = interpol(nhis,narr,n_absc) ; Interpolate onto common scale
     phys[ii].npdf /= total(phys[ii].npdf) ; Normalize to unit integral prob.
     cgPlot,n_absc,phys[ii].npdf,/xst
     vline,4,color='blue'
     vline,alog10(750),color='blue'
     vline,/h,max(phys[ii].npdf)/4.d,color='dark orchid'
     IF NJ GE 2 THEN message,'Fire 4',/inf
     
     plothist,obj.n20k,bin=median(obj.n20k)/1.d3,narr,nhis,/noplot
     phys[ii].n20k = omni_compute_ebars(nhis,narr,/logcut,NJ=nj)
     IF NJ GE 2 THEN message,'Fire 5',/inf
     
     
     
     ;;&&&&&&&& Print Line for LaTeX Table &&&&&&&&
     lprad = 10.^phys[ii].prad[0] * $
             abs([1., (1. - 10.^(-phys[ii].prad[1])),$
                  (1. - 10.^(-phys[ii].prad[2]))])
     
     ;; print,phys[ii].prad
     ;; print,lprad
     ;; STOP
     
     
     prad_nden = ~unresolved ? $
                 string( lprad[0] GE 0. ? '' : '',$
                         lprad[0],lprad[1],lprad[2],$
                         phys[ii].nden[0],phys[ii].nden[1],phys[ii].nden[2],$
                         format="(A0,'$',F0.2,'_{-',F0.2,'}^{+',F0.2,"+$
                         "'}$ & $',F0.2,'_{-',F0.2,'}^{+',F0.2,'}$')") : $
                 '\nodata & \nodata'

     catstr = string(s[ii].cnum,s[ii].glon,$
                       s[ii].glat GE 0 ? '\phs' : '',s[ii].glat,$
                       s[ii].(iflux),s[ii].(ieflux),thetar,$
                       pvec[ii].stat.duse[0]/1.d3 LT 10. ? '\phn' : '',$
                       pvec[ii].stat.duse[0]/1.d3,$ 
                       pvec[ii].stat.duse[1]/1.d3,pvec[ii].stat.duse[2]/1.d3,$
                       phys[ii].mass[0],phys[ii].mass[1],phys[ii].mass[2],$
                       prad_nden,$
                       format="(I0,' & ',F0.3,' & ',A0,'$',F0.3,'$ & $',F0.2,"+$
                       "'(',F0.2,')$ & ',A0,' & ',A0,'$',F0.2,"+$
                       "'_{-',F0.2,'}^{+',F0.2,'}$ & $',F0.2,'_{-',F0.2,"+$
                       "'}^{+',F0.2,'}$ & ',A0,' \\')")
     

     printf,lun,catstr
     IF ii LT 10 THEN BEGIN
        printf,lun2,catstr
        flush,lun2              ; Force writing to file for the impatient user
     ENDIF
     flush,lun                  ; Force writing to file for the impatient user
     
     ;; Print ETD statistics for the impatient user
     etd = start_t + (systime(1) - start_t) / (double(ii+1)/nkp)
     message,'ETD:  '+systime(0,etd),/inf
     
  ENDFOR                        ; End of loop over sources
  
  close,lun
  free_lun,lun
  close,lun2
  free_lun,lun2
  
  save,phys,mmin,mmax,rmin,rmax,nmin,nmax,$
       filename='./masses/save_files/physics.sav',/ver
  
  ;;&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
  ;; Make plots!
  
  myps,'./masses/plots/maxvals.eps',xsize=14,ysize=10
  multiplot_xm,[3,2],gap=0.035,mpcharsize=1.0,/doxaxis
  
  plothist,charsize=1.0,mmax,xtit='Mass'
  
  multiplot,/doyaxis,/doxaxis
  
  plothist,charsize=1.0,rmax,xtit='Radius'
  
  multiplot,/doyaxis,/doxaxis
  
  plothist,charsize=1.0,nmax,xtit='Number Density'
  
  multiplot,/doyaxis,/doxaxis
  
  plothist,charsize=1.0,mmin,xtit='Mass'
  
  multiplot,/doyaxis,/doxaxis
  
  plothist,charsize=1.0,rmin,xtit='Radius'
  
  multiplot,/doyaxis,/doxaxis
  
  plothist,charsize=1.0,nmin,xtit='Number Density'
  
  myps,/done,/mp
  
END
